This is an analysis of the PTSD model “complete” simulations (“Simulation 3”).

# Simulate ACT-R

decay <- function(time, t0=0, rate=.5) {
    if (time <= t0) {
      NA
    } else {
      (time - t0)**(-rate)
    }
}

vdecay <- Vectorize(decay)

time <- seq(1, 100, 1)

t1 = vdecay(time, t0=1)
t2 = vdecay(time, t0=20)
t3 = vdecay(time, t0=55)
t4 = vdecay(time, t0=65)


df <-data.frame(Time=time, 
           trace1 = vdecay(time, t0=1),
           trace2 = vdecay(time, t0=20),
           trace3 = vdecay(time, t0=55),
           trace4 = vdecay(time, t0=65))


ldf <- df %>%
  # mutate(Total = if_else(is.na(trace1), 0, trace1) +
  #          if_else(is.na(trace2), 0, trace2) +
  #          if_else(is.na(trace3), 0, trace3) +
  #          if_else(is.na(trace4), 0, trace4))  %>%
  pivot_longer(cols = c("trace1", "trace2", "trace3", "trace4"),
               names_to = "Trace", values_to="Activation")

ldf$Trace <- fct_recode(ldf$Trace, 
                        "Trace 1" = "trace1",
                        "Trace 2" = "trace2",
                        "Trace 3" = "trace3",
                        "Trace 4" = "trace4")

ldf <- ldf %>%
  add_column(Source="Individual Traces")

t1[is.na(t1)] <-0
t2[is.na(t2)] <-0
t3[is.na(t3)] <-0
t4[is.na(t4)] <-0
total <- log(t1+t2+t3+t4)

totaldf <- tibble(Time=time, 
                  Activation = total,
                  Trace = "Memory",
                  Source="Memory")

ldf <- rbind(ldf, totaldf)

ldf$Source <- factor(ldf$Source, levels = c("Memory", "Individual Traces"))

ggplot(ldf, 
       aes(x=Time, y=Activation, col=Trace)) +
  geom_line(size=1,
            alpha=1) +
            #linetype="dashed") +
  # stat_summary(geom="line", fun = sum, 
  #              col="black", 
  #              alpha=0.5,
  #              position = position_dodge()) +
  # #scale_y_log10() +
  facet_wrap(~ Source, ncol=1,
             scales = "free_y") +
  #ylim(0, 1.6) +
  #ylab(expression(paste("Availability of Memory ", italic(m))))+
  ylab("Activation") +
  scale_color_viridis("", option="plasma", discrete=T, end = .8, direction=-1) +
  theme_pander()

## Warning: Removed 141 row(s) containing missing values (geom_path).

The distribution of Gamma events is this:

Frequency: 2.1, 175 (shape, scale)
Rumination, 6, 100

t <- 0:(24*60)
E <- dgamma(t, shape=2.1, rate=1/175)
R <- dgamma(t, shape=6, rate=1/100)

dists <- data.frame(Time=t, Events = E, Rumination = R) %>%
  pivot_longer(names_to = "Type", 
               values_to = "Probability",
               cols=c("Events", "Rumination"))

ggplot(dists, aes(x=Time, y=Probability, col=Type, fill=Type)) +
  geom_line() +
  geom_ribbon(data=dists, aes(x=Time, ymax=Probability, ymin=0), alpha=0.25) +
  scale_x_continuous(breaks=60*c(0, 4, 8, 12, 16, 20),
                     limits=c(0, 24*60),
                     labels=c("8:00am", 
                              "12:00pm",
                              "4:00pm",
                              "8:00pm",
                              "12:00am",
                              "4:00am")) +
  scale_color_aaas() +
  scale_fill_aaas() +
  ggtitle("Probability of Retrievals During the Day") +
  theme_pander()

Load and Transform the Data

First, we are going to load the “aggregated” version of the simulations, with the events already aggregated by day.

CREATE_AGGREGAGATED = F

if (CREATE_AGGREGAGATED) {
  d <- read_csv("../simulations/simulations3.csv")
  #d <- read_csv("simulations3.csv", col_types = cols())
  d <- d %>%
    mutate(Day = floor(Time / (60*60*24)) - 100,
           Hour = floor((d$Time / 3600) %% 24),
           RuminationFrequency=as_factor(RuminationFrequency))
  names(d)[18] <- "W"
  
  d <- d %>% 
    dplyr::select(Run, Day, PTEV, Gamma, PTES, W, 
                  NumAttributes, RuminationFrequency, 
                  MemoryEntropy, ChunkSimilarity, Traumatic, ChunkV) %>%
    filter(Day > -20)
  
  a <- d  %>%
    group_by(Run, Day, PTEV, PTES, Gamma, 
             NumAttributes, RuminationFrequency, W) %>%
    summarise(MemoryEntropy = mean(MemoryEntropy),
              ChunkSimilarity = mean(ChunkSimilarity),
              PTraumatic = mean(Traumatic),
              CTraumatic = sum(Traumatic),
              ChunkV = mean(ChunkV))
  
  write_csv(x=a, file = "../simulations/simulations3_aggregated2.csv")
}

a <- read_csv(gzfile("../simulations/simulations3_aggregated2.csv.gz"),
              col_types = cols())

Then, we are going to rename the simulation variables to make them consistent with the names used in published papers:

a <- a %>%
  rename(I = PTEV) %>%
  rename(gamma = Gamma) %>%
  rename(A = NumAttributes) %>%
  rename(C = PTES) %>%
  rename(R = RuminationFrequency)

Recovery Trajectories

To identify a recovery trajectory, we need to first define a classification function. Here, we are going to use the following algorithm:

First the model calculates three average values:

The mean probability P(baseline) of retrieving intrusive memories in the 10 days preceding the PTE, or baseline period. By definition, this is always zero.
The mean probability P(acute) of retrieving intrusive memories in the 10 days following the PTE, or during the acute period.
The mean probability P(chronic) of retrieving intrusive memories in the last ten days of the second month after the PTE (i.e, days 51-60), or the chronic period.

Then, the model calculates two statistical tests:

A t-test between the baseline and acute period, or acute test; and
A t-test between the acute-period and the chronic period, or chronic test.

And finally we apply the following decision tree:

If the acute test is significant at p < .05 (Pacute > Pbaseline), then

1.1. If the chronic test is also significant at p < .05, and Pchronic > Pacute, we classify the trajectory as Delayed.

1.2. If the chronic test is significant at p < .05 but P(chronic) < P(acute), then we classify the trajectory as Recovery;

1.3. If the chronic test is non-significant, then P(chronic) ~ P(acute) , and we classify the trajectory as Chronic.
If, instead, the acute test is not significant and P(acute) ~ P(baseline), then:

2.1. If the chronic test is significant at p < .05, then P(chronic) > P(acute), and we classify the trajectory as Delayed.

2.2 If the chronic test is also not significant, then P(chronic) ~ P(baseline) ~ P(acute), and we classify the trajectory as Resilient.

# Classifies people based on trajectory:
# Chronic = 1
# Recovery =2
# Delay    =3
# resilient= 4
#
classify <- function(days) {
  g1 <- days[10:19]
  g2 <- days[21:30]
  g3 <- days[70:79]
  t12 <- t.test(g1, g2)
  t23 <- t.test(g2, g3)
  category <- 0
  if (!is.na(t12$p.value) & t12$p.value < 0.05) {
    # t1 < t2
    if (!is.na(t23$p.value) & t23$p.value < 0.05) {
      # t2 <> t3

      if (!is.na(t23$statistic) & t23$statistic > 0) {
        # t1 < t2, t2 > 3
        category <- "Recovery" #2  # Recovery
      } else {
        # t1 < t2, t2 < t3
        category <- "Delayed" # 3 # Delayed
      }
    } else {
      # t1 < t2, t2 = t3
      category <- "Chronic" # 1 # Chronic
    }
  } else {
    # t1 == t2
    if (!is.na(t23$p.value) & t23$p.value < 0.05) {
      # t1 == t2, t2 < t3
      category <- "Delayed" #3 # Delayed
    } else {
      # t1 == t2, t2 = t3
      category <- "Resilient" # 4 # Resilient
    }
  }
  category
}

patterns <- a %>%
  group_by(Run, I, C, A, R, W, gamma) %>%
  dplyr::summarize(Trajectory = classify(PTraumatic))

## `summarise()` has grouped output by 'Run', 'I', 'C', 'A', 'R', 'W'. You can override using the `.groups` argument.

patterns <- patterns %>% ungroup

Finally, let’s assign each data point in the main averages with the corresponding trajectory, and re-order the trajectories in order of severity:

a <- inner_join(a, patterns)

## Joining, by = c("Run", "I", "C", "gamma", "A", "R", "W")

a$Trajectory <- factor(a$Trajectory,
                      levels = c("Recovery", "Delayed", "Resilient",  "Chronic"))

Chronic Probability of Traumatic Intrusions

The other “behavioral” measure is the incidence of traumatic memory intrusions. As an aggregate measure, we will consider only the averages in the “Chronic” period, which we will be taking as indicative of long-term consequences of the PTSD.

traumatic <- a %>%
  filter(Day > 1) %>%
  group_by(Run, I, C, R, W, gamma) %>%
  summarize(PTraumatic = mean(PTraumatic),
            CTraumatic = sum(CTraumatic))

## `summarise()` has grouped output by 'Run', 'I', 'C', 'R', 'W'. You can override using the `.groups` argument.

At this point, we can visualize the mean aggregated timecourse of intrusions for each trajectory.

K = rev(brewer.pal(5, "Purples"))

TOTAL <- 28800  # Param combos

traj_total <- a %>%
  filter(Day == 30) %>%
  group_by(Trajectory) %>%
  summarize(N = length(CTraumatic),
            Prop = length(CTraumatic)/TOTAL,
            MeanC = mean(CTraumatic))



ggplot(filter(a, I != 1), 
       aes(x=Day, y=CTraumatic, col=Trajectory, fill=Trajectory)) +
  ggtitle("Probability of Memory Intrusion by Trajectory") +
  #geom_vline(data=props, aes(xintercept=-0.5)) +
  annotate("segment", x=-0, xend=-0,
           y=-Inf, yend=Inf, col="black", lty=2, size=1) +
  annotate("rect", xmin=50, xmax=60,
           ymin=-Inf, ymax=Inf, fill=K[1], alpha=0.3)+
  annotate("rect", xmin=1, xmax=10,
           ymin=-Inf, ymax=Inf, fill=K[2], alpha=0.3)+
  annotate("rect", xmin=-10, xmax=-1,
           ymin=-Inf, ymax=Inf, fill=K[3], alpha=0.3)+
  
  annotate("text", x= -5, y=15, label="Base")+
  annotate("text", x= 5, y=15, label="Acute")+
  annotate("text", x= 55, y=15, label="Chronic")+
  
  geom_text(data=traj_total, aes(x= 30, y = MeanC, 
                                 label=percent(Prop, accuracy = 0.1)),
            vjust=-.5, show.legend=F) +
  
  stat_summary(fun.data = mean_se, geom="line") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="ribbon", alpha = 0.25, col=NA) +
  coord_cartesian(ylim=c(0, 15)) +
  ylab("Total Number of Memory Intrusions") +
  scale_color_d3() +
  scale_fill_d3() +
  theme_pander()

The results indicate that our classification was correct; the trajectories behave as planned.

Finding a baseline to compare the model to existing literature

Since we have no idea how typical the model parameters are of the existing patients and their cases, we need to provide some sort of baselines. To do so, we will use the proportion of trajectories that were given by Galatzer-Levy and Bonanno in their meta-analysis. They are:

Resilient, 0.657
Recovery, 0.208
Chronic, 0.106
Delayed, 0.089

Note that they do not normalize to one, due to the lack of some trajectories in some studies. To use them, we will need to make sure to normalized them first, which can be done at time of testing.

T_observed <- c(0.657, 0.208, 0.106, 0.089)
names(T_observed) <- c("Resilient", "Recovery", "Chronic", "Delayed")

We can now plot the percentages and the patterns by trajectory:

props <- patterns %>%
  group_by(Trajectory, 
           I, 
           gamma, 
           C, 
           W,
           R,
           #A,
           ) %>%
  dplyr::summarize(N = length(Trajectory)) %>%
  dplyr::mutate(Prop = N/50)

## `summarise()` has grouped output by 'Trajectory', 'I', 'gamma', 'C', 'W'. You can override using the `.groups` argument.

props %>%
  pivot_wider(id_cols=c(I, 
                        gamma, 
                        C,
                        W,
                        R,
                        #A,
                        ), 
              names_from = Trajectory,
              values_from = N,
              values_fill=0) -> test

mychi <- function(vals) {
  chisq.test(vals, 
             p=T_observed, 
             rescale.p=T)$statistic
}

mycorr <- function(vals) {
  cor(vals, T_observed)
}

myrmse <- function(vals) {
  v <- 50 * (T_observed/sum(T_observed))
  sqrt(mean((vals - v)^2))
}


lprops <- test %>%
  pivot_longer(cols=c("Resilient", "Recovery", "Chronic", "Delayed"),
               names_to = "Trajectory",
               values_to = "N")

lprops$Trajectory <- factor(lprops$Trajectory,
                      levels = c("Recovery", "Delayed", "Resilient",  "Chronic"))

# Suppresses X-squared approximation warnings
old.warn <- options()$warn
options(warn = -1)

atest <- lprops %>%
  group_by(I, 
           gamma, 
           C, 
           W,
           R,
           #A,
           ) %>%
  summarize(Chi = mychi(N), r=mycorr(N), RMSE=myrmse(N))

## `summarise()` has grouped output by 'I', 'gamma', 'C', 'W'. You can override using the `.groups` argument.

# resets warnings
options(warn = old.warn)

test <- inner_join(test, atest)

## Joining, by = c("I", "gamma", "C", "W", "R")

Now, we can identfy the baseline as the condition with the minimum \(\chi^2\)-squared test (or the minimum RMSE, they are the same)

This shows how correlation coefficient and \(\chi^2\) are related, while the RMSE does not really capture either of them:

ggplot(test, aes(x=Chi, y=r, col=RMSE)) + 
  geom_point(size=4, alpha=0.5) + 
  scale_color_viridis(option="plasma")+
  ylab(expression(italic(r))) +
  xlab(expression(chi^2)) +
  ggtitle("Relationship between Error Measures in Trajectories") +
  theme_pander()

Now, let’s identify a “baseline” condition—the one that most resembles the distribution of trajectories identified by Galatzer-Levy.

baseline <- test %>%
  filter(Chi == min(test$Chi))

baseline %>%
  kable(digits=3) %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

I	gamma	C	W	R	Chronic	Delayed	Recovery	Resilient	Chi	r	RMSE
40	0.9	0.25	8	0	15	3	15	67	7.464	0.985	18.875

And here is a visual representation of the four trajectories at that baseline (smoothed, because there is too much day-to-day variance in this small sample).

baseline_a <- a %>%
  filter(I == baseline$I,
         W == baseline$W,
         C == baseline$C,
         R == baseline$R,
         gamma == baseline$gamma,
         )

base_total <- baseline_a %>%
  filter(Day == 15) %>%
  group_by(Trajectory) %>%
  summarize(N = length(CTraumatic),
            Prop = length(CTraumatic)/100,
            MeanC = mean(CTraumatic))

         
ggplot(baseline_a, 
       aes(x=Day, y=CTraumatic, col=Trajectory, fill=Trajectory)) +
  geom_smooth(method="loess", span=0.2) +
  ggtitle("Trajectories for Best-Fitting Parameters") +
  scale_color_d3() +
  geom_text(data=base_total, aes(x= c(25, 10, 20, 10), y = MeanC, 
                                 label=percent(Prop, accuracy = 0.1)),
            vjust=-.5, show.legend=F) +
  ylab("Total Number of Memory Intrusions") +
  geom_vline(data=props, aes(xintercept=-0.25)) +
  scale_fill_d3() +
  theme_pander()

## `geom_smooth()` using formula 'y ~ x'

The Chonic and Delayed trajectories, being the least common, exhibit significant ups and downs, hence the need for Loess smoothing. Nonetheless, the trajectories remain visible and clearly identifiable.

But how good is the baseline condition, in terms of being closed to Galatzer-Levy’s estimated pooled values? We can just examine the significance of the \(\chi^2\) test.

baseline_vals <- baseline[c("Resilient", "Recovery", "Chronic", "Delayed")]
chi <- chisq.test(baseline_vals, p=T_observed, rescale.p = T)

The result is not too bad—A bit dissimilar, but not too far either, and not statistically significant (\(\chi^2\) = 7.463523, p = 0.058503)).

The main difference is that the delayed onset trajectory has a higher proportion than in the review. However, Galatzer-Levy and Bonanno themselves note that the resilient (M= 0.66) delayed onset trajectories (M = 0.099) has a higher prevalence in studies of single-point event. If these corrected prevalences are taken into account, the chi squared reduces further.

Factors affecting the distribution of trajectories

Although theoretically possible, a full five-way analysis of model’s factors and interactions would be statistically uninformative. Instead, we will examine the model’s outcomes and behaviors with the respect to the known main effects reported in the literature.

We can now examine, one by one, how the different factor in the model affect the distribution of trajectories.

To examine these factors, we need to use logistic regression. And, to do so, we need to create a number of dummy variables specifying a binary status for each trajectory in patterns.

patterns <- patterns %>%
  mutate(Resilient = if_else(Trajectory == "Resilient", 1, 0),
         Chronic = if_else(Trajectory == "Chronic", 1, 0),
         Recovery = if_else(Trajectory == "Recovery", 1, 0),
         Delayed = if_else(Trajectory == "Delayed", 1, 0))

External factors: PTE Intensity and Context Similarity

There are two main external factors in the model: The intensity of the PTE (I) and the contextual similarity (C). They will be both examined here.

Effects on Timecourses

First, let’s load additional values of I and C in addition to the ones in the baseline.

IC_a <- a %>%
  filter(#I == baseline$I,
    I != 1,
         W == baseline$W,
         #C == baseline$C,
         R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )

IC_patterns <- patterns %>%
    filter(#I == baseline$I,
    I != 1,
         W == baseline$W,
         #C == baseline$C,
         R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )

IC_trajectories <- lprops %>%
    filter(#I == baseline$I,
    I != 1,
         W == baseline$W,
         #C == baseline$C,
         R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )


IC_traumatic <- traumatic %>%
      filter(#I == baseline$I,
         I != 1,
         W == as.numeric(baseline$W),
         #C == baseline$C,
         R == as.numeric(baseline$R),
         gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

First, let’s look at the time course of intrusions:

ggplot(IC_a, 
       aes(x=Day, y=CTraumatic, col=as.factor(I), fill=as.factor(I))) +
  #geom_smooth(method="loess", span=0.2) +
  stat_summary(fun.data = mean_se, geom="line") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="ribbon", alpha = 0.25, col=NA) +
  
  facet_wrap(~C, labeller = label_bquote(italic(C) == .(C))) + #label=label_both) +
  labs(col=expression(italic(I)[PTE]), 
       fill=expression(italic(I)[PTE]), 
       facet="Z") +
  #coord_cartesian(ylim=c(0, 0.5)) +
  ggtitle("Effects of PTE Intensity") +
  ylab("Number of Traumatic Memory Intrusions") +
  
  scale_color_viridis(discrete=T, option="plasma", end=0.8) +
  scale_fill_viridis(option = "plasma", discrete = T, end=0.8) +
  geom_vline(data=props, aes(xintercept=-0.15)) +
  theme_pander()

The effects of I and C on the probabilities of recall during the Chronic period can be understood with a linear model analysis

IC_mod <- glm(CTraumatic ~ I * C, family="poisson", IC_traumatic)
summary(IC_mod) %>%
  xtable %>%
  kable %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Estimate	Std. Error	z value
(Intercept)	-0.7993158	0.0439928	-18.16923
I	0.0842682	0.0008007	105.23676
C	7.8729269	0.0636742	123.64387
I:C	-0.0739858	0.0011709	-63.18745

Effects on Trajectories

And here are the changes in trajectories.

ggplot(data=IC_trajectories, aes(x="", y=N/100, fill=Trajectory)) +
  geom_bar(stat = "identity", col="grey", width=1) +
  #facet_grid(C~I, labeller = label_both) +
  facet_grid(C~I, labeller = label_bquote(
    rows = italic(C) == .(C),
    cols = italic(I)[PTE] == .(I))) +
  coord_polar("y", start=0) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 2L)) +
  scale_fill_d3() +
  scale_color_d3() +
  ylab("PTE Intensity") +
  xlab("Contextual Similarity") +
  ggtitle("Effect of External Factors on Trajectories") +
  geom_text_repel(aes(label=percent(N/100, .1)), col="white",
            position=position_stack(vjust=0.5 ), size=3)+
  theme_pander() +
  theme(axis.ticks = element_blank(), 
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        legend.position = "bottom")

Are the change in trajectories meaningful? To do so, we first examine the trajectories using a MANOVA:

man <- manova(cbind(Resilient, Recovery, Delayed) ~ I * C, IC_patterns)
tidy(man) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

term	df	pillai	statistic	num.df	den.df	p.value
I	1	0.2505389	133.0483	3	1194	0
C	1	0.5364246	460.5442	3	1194	0
I:C	1	0.0589542	24.9337	3	1194	0
Residuals	1196	NA	NA	NA	NA	NA

Then, we examine each category independently, we will perform a logistic regression on dummy binary variables that categorize each run of the model as belonging to that category, or not.

In the case of the Resilient and Recovery trajectories, both intensity I, Context C are significant, but their interaction is not. This implies that the effects of I and C are additive.

Furthermore, for the Resilient trajectories, the coefficients are negative, implying that the proportion of Recovery and Resilient trajectories decreases as the intensity and the context similarity increase.

Resilient Trajectory

mylogit <- glm(Resilient ~ (I*C), family="binomial", data=IC_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	8.0726911	0.7396808	10.9137497	0.0000000
I	-0.1042898	0.0141387	-7.3761937	0.0000000
C	-10.3446128	1.4043857	-7.3659343	0.0000000
I:C	-0.0144672	0.0322948	-0.4479746	0.6541716

Recovery Trajectory

For Recovery trajectories, however, the coefficients are significantly positive, implying that the number of cases increase with the severity of intensity and similarity. This is likely due to a greater number of resilient models developing intrusive memories following PTE under greater contextual similarity.

mylogit <- glm(Recovery ~ (I*C), family="binomial", data=IC_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-6.5266334	0.7852632	-8.3113963	0.0000000
I	0.0757766	0.0147198	5.1479313	0.0000003
C	5.0969156	1.2604091	4.0438580	0.0000526
I:C	-0.0085358	0.0247694	-0.3446121	0.7303860

Chronic Trajector

In the case of Chronic and Delayed trajectories, however, the pattern is different. In these cases, both the main factors and their interaction are significant. Furthermore, in both cases, the main factors have positive coefficients but their interaction is negative. This suggests that the effects of one factor decrease the size of the other factor increases.

mylogit <- glm(Chronic ~ (I*C), family="binomial", data=IC_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value
(Intercept)	-6.2050051	0.6318950	-9.819677
I	0.0760648	0.0121663	6.252087
C	9.4454215	1.0400816	9.081424
I:C	-0.1332185	0.0207257	-6.427706

Delayed Trajectory

mylogit <- glm(Delayed ~ (I*C), family="binomial", data=IC_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-9.2380905	1.7014060	-5.429680	0.0000001
I	0.1004351	0.0311299	3.226324	0.0012539
C	11.4814812	2.4891026	4.612699	0.0000040
I:C	-0.2098212	0.0489880	-4.283111	0.0000184

Internal Factors: Working Memory and Recollection Vividness

The internal factors are two, differences in WMC and differences in vividness.

This can be related to the difference in children.

Wg_a <- a %>%
  filter(I == baseline$I,
         #W == baseline$W,
         C == baseline$C,
         R == baseline$R,
         #gamma == baseline$gamma,
         #A == baseline$A
         )

Wg_test <- lprops %>%
    filter(I == baseline$I,
         #W == baseline$W,
         C == baseline$C,
         R == baseline$R,
         #gamma == baseline$gamma,
         #A == baseline$A
         )

Wg_patterns <- patterns %>%
    filter(I == baseline$I,
         #I != 1,
         #W == as.numeric(baseline$W),
         C == baseline$C,
         R == as.numeric(baseline$R),
         #gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

Wg_traumatic <- traumatic %>%
    filter(I == baseline$I,
         #I != 1,
         #W == as.numeric(baseline$W),
         C == baseline$C,
         R == as.numeric(baseline$R),
         #gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

Effects on the Timecourse of Intrusive Memories

Here we visualize the change in timecourses.

Wg_a$gamma <- as_factor(Wg_a$gamma)
ggplot(Wg_a, 
       aes(x=Day, y=CTraumatic, col=gamma, fill=gamma)) +
#  geom_point(alpha=0.1) +
  stat_summary(fun.data = mean_se, geom="line") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="ribbon", alpha = 0.25, col=NA) +
  
  #geom_smooth(method="loess", span=0.2) +
  facet_wrap(~W, ncol=3,
             labeller  = label_bquote(
               rows = italic(W) == .(W))) +
  geom_vline(data=props, aes(xintercept=-0.15)) +
  ylab("Number of\nTraumatic Memory Intrusions") +
  #coord_cartesian(ylim=c(0, 0.75)) +
  ggtitle("Effects of Idiographic Parameters") +
  scale_color_viridis(expression(paste(gamma, " =")), discrete=T, 
                      option="plasma", begin = 0, end=0.9) +
  scale_fill_viridis(expression(paste(gamma, " =")), option = "plasma", 
                     discrete = T, begin = 0, end=0.9) +
  theme_pander() +
  theme(legend.position = "bottom")

It is striking to notice that, within baseline levels of I and C, the effects of working memory are incredibly large. This is likely due to the fact that the parameters were not properly calibrated.

And here is the corresponding analysis. Like before, the main factors W and \(\gamma\) are significant. However, this time, their interaction is not, suggesting that their effects are additive.

Wg_mod <- glm(CTraumatic ~ W * gamma, family="poisson", Wg_traumatic)
summary(Wg_mod) %>%
  xtable %>%
  kable %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-4.6088943	0.2487333	-18.5294648	0.0000000
W	-0.5902732	0.0547969	-10.7720158	0.0000000
gamma	15.0052321	0.2687441	55.8346535	0.0000000
W:gamma	-0.0295063	0.0592083	-0.4983467	0.6182397

Effects on Trajectories

And here are the changes in trajectories.

ggplot(data=Wg_test, aes(x="", y=N/100, fill=Trajectory)) +
  geom_bar(stat = "identity", col="white", width=1, alpha=1) +
  facet_grid(gamma~W, 
             labeller  = label_bquote(
               rows = gamma == .(gamma),
               cols = italic(W) == .(W))) +
  coord_polar("y", start=0) +
  #scale_fill_d3() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 2L)) +
  scale_fill_d3() +
  xlab("Recollection Vividness") +
  ylab("Cognitive Control") +
  ggtitle("Effect of Idiographic Factors on Trajectories") +
  geom_text_repel(aes(label=percent(N/100, .1)),
            position=position_stack(vjust=0.5), size=3, col="white") +
  #geom_text_repel() +
  theme_pander() +
  theme(axis.ticks = element_blank(), 
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        legend.position = "bottom")

As before, we can start with a MANOVA:

man <- manova(cbind(Resilient, Recovery, Delayed) ~ W * gamma, Wg_patterns)
tidy(man) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

term	df	pillai	statistic	num.df	den.df	p.value
W	1	0.6151552	476.33812	3	894	0
gamma	1	0.0818391	26.56184	3	894	0
W:gamma	1	0.0478437	14.97384	3	894	0
Residuals	896	NA	NA	NA	NA	NA

The MANOVA shows that both factors and their interactions are significant. We can now carry out a trajectory by trajectory analysis.

First, we can examine the effects on the Resilient trajectory. Because of the massive amount of positive cases when W = 12, we are going to use the bigglm function, which handles cases better.

Resilient Trajectory

mylogit <- bigglm(Resilient ~ W * gamma, 
                  family = binomial("logit"),
                  data = Wg_patterns)

## Warning in bigglm.function(formula = formula, data = datafun, ...): ran out of
## iterations and failed to converge

summary(mylogit)

## Large data regression model: bigglm(Resilient ~ W * gamma, family = binomial("logit"), data = Wg_patterns)
## Sample size =  900 
## failed to converge after 8  iterations
##                  Coef      (95%      CI)      SE      p
## (Intercept)   88.4971    8.8649 168.1293 39.8161 0.0262
## W             -9.8782  -19.8616   0.1052  4.9917 0.0478
## gamma       -116.3745 -215.7721 -16.9768 49.6988 0.0192
## W:gamma       13.3382    0.8849  25.7914  6.2266 0.0322

Surprisingly, there is no significant effect of either parameter on the Resilient trajectory. This is likely due to the massive amount of variance, with the Resilient Trajectory going from being almost absent when W = 4 to being the only possible trajectory when W = 12.

Chronic Trajectory

We can then examine the Chronic trajectory:

mylogit <- glm(Chronic ~ (W * gamma), family=binomial("logit"), 
               data=Wg_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-11.1345718	3.7901285	-2.9377822	0.0033057
W	0.4044620	0.6134173	0.6593586	0.5096655
gamma	13.9799411	4.2315576	3.3037341	0.0009541
W:gamma	-0.9254449	0.6835085	-1.3539625	0.1757483

Recovery Trajectory

Then, we on the Recovery trajectory.

mylogit <- glm(Recovery ~ (W * gamma), family=binomial("logit"), 
               data=filter(Wg_patterns, W<13))
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	14.111736	4.3608454	3.236009	0.0012121
W	-2.075010	0.8391947	-2.472620	0.0134127
gamma	-13.420202	4.9093783	-2.733585	0.0062649
W:gamma	1.638059	0.9425624	1.737878	0.0822323

Delayed Trajectory

The Delayed trajectory:

mylogit <- glm(Delayed ~ (W * gamma), 
               family=binomial("logit"), 
               data=filter(Wg_patterns, W < 20) )
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-1.481706	8.516561	-0.1739794	0.8618817
W	-2.680543	1.722150	-1.5565091	0.1195871
gamma	1.304849	9.171375	0.1422741	0.8868635
W:gamma	2.549777	1.841799	1.3843945	0.1662377

Comorbidty: Effects of Rumination

R_a <- a %>%
  filter(I == baseline$I,
         W == baseline$W,
         C == baseline$C,
         #R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )

R_test <- lprops %>%
    filter(I == baseline$I,
         W == baseline$W,
         C == baseline$C,
         #R == baseline$R,
         gamma == baseline$gamma,
         #A == baseline$A
         )

R_patterns <- patterns %>%
    filter(I == baseline$I,
         #I != 1,
         W == as.numeric(baseline$W),
         C == baseline$C,
         #R == as.numeric(baseline$R),
         gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

R_traumatic <- traumatic %>%
    filter(I == baseline$I,
         #I != 1,
         W == as.numeric(baseline$W),
         C == baseline$C,
         #R == as.numeric(baseline$R),
         gamma == as.numeric(baseline$gamma),
         #A == baseline$A
         )

Effects on the Timecourse of Intrusive Memories

Here we visualize the change in timecourses.

ggplot(R_a, 
       aes(x=Day, y=CTraumatic, col=as.factor(R), fill=as.factor(R))) +
  stat_summary(fun.data = mean_se, geom="line") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="ribbon", alpha = 0.25, col=NA) +
  geom_vline(data=props, aes(xintercept=-0.15)) +

  #geom_smooth(method="loess", span=0.2) +
  #facet_wrap(~R) +
  #coord_cartesian(ylim=c(0, 0.75)) +
  ylab("Number of\nTraumatic Memory Intrusions") +
  #coord_cartesian(ylim=c(0, 0.75)) +
  ggtitle("Effects of Rumination") +
  scale_color_viridis(expression(paste(italic(R), " =")), discrete=T, option="plasma", end=0.8) +
  scale_fill_viridis(expression(paste(italic(R), " =")), option = "plasma", discrete = T, end=0.8) +
  theme_pander() +
  theme(legend.position = "bottom")

The effects of Rumination are huge. This is likely due to the fact that the parameters were not properly calibrated.

And here is the corresponding analysis. Like before, the main factors W and \(\gamma\) are significant. However, this time, their interaction is not, suggesting that their effects are additive.

R_mod <- glm(CTraumatic ~ R, family="poisson", R_traumatic)
summary(R_mod) %>%
  xtable %>%
  kable %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	4.0167433	0.0189797	211.6333	0
R	0.1350059	0.0009804	137.7115	0

Effects on Trajectories

And here are the changes in trajectories.

ggplot(data=R_test, aes(x="", y=N/100, fill=Trajectory)) +
  geom_bar(stat = "identity", col="white", width=1, alpha=1) +
  facet_wrap(~ R, labeller = label_bquote(
               rows = italic(R) == .(R))) +
  coord_polar("y", start=0) +
  #scale_fill_d3() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 2L)) +
  scale_fill_d3() +
  ylab("Rumination") +
  xlab("") +
  ggtitle("Effect of Rumination on Trajectories") +
  geom_text_repel(aes(label=percent(N/100, .1)),
            position=position_stack(vjust=0.5), size=3, col="white") +
  #geom_text_repel() +
  theme_pander() +
  theme(axis.ticks = element_blank(), 
        axis.text.y = element_blank(),
        axis.text.x = element_blank())

As before, we can start with a MANOVA:

man <- manova(cbind(Resilient, Recovery, Delayed) ~ R, R_patterns)
tidy(man) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

term	df	pillai	statistic	num.df	den.df	p.value
R	1	0.4574771	55.09169	3	196	0
Residuals	198	NA	NA	NA	NA	NA

The MANOVA shows that rumination is significant. We can now carry out a trajectory by trajectory analysis.

First, we can examine the effects on the Resilient trajectory.

mylogit <- glm(Resilient ~ R, 
               family = "binomial",
               data = R_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	0.7081851	0.2126697	3.329976	0.0008685
R	-0.2092142	0.0311797	-6.709947	0.0000000

We can then examine the Chronic trajectory:

mylogit <- glm(Chronic ~ R, family="binomial", data=R_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-1.7346011	0.2800560	-6.193765	0.0000000
R	0.0622526	0.0173836	3.581108	0.0003421

The, we on the Recovery trajectory.

mylogit <- glm(Recovery ~ R, family="binomial", data = R_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-1.7346010	0.2800209	-6.194541	0
R	0.0967636	0.0172348	5.614430	0

And, finally, The Delayed trajectory:

mylogit <- glm(Delayed ~ R, 
               family="binomial", 
               data=R_patterns)
summary(mylogit) %>%
  xtable() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("hover", "striped"))

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-3.4760987	0.5862093	-5.9297915	0.0000000
R	0.0149022	0.0388606	0.3834797	0.7013641

Neurobiology of PTSD

Finally, we will examine how the model captures some known effects of the neurobiology of PTSD and, in particular, the size of the hippocampus.

Hippocampus Size

First, we define a baseline condition: No PTE (I = 1) and \(\gamma\) = 0.8.

H1 <- a %>%
  filter(I == 1 & gamma == 0.9 & C==0.25 & R == 0 & Day > 49) %>%
  dplyr::select(MemoryEntropy) %>%
  summarise(H=mean(MemoryEntropy)) %>%
  as.double()

Then, we define a function that calculates memory entropy on the last N days of a timeseries.

h_entropy <- function(timeseries, N = 10) {
  L <- length(timeseries)
  mean(timeseries[L-N:N])
}

With this in place, we can now calculate the predicted decrease in hippocampus size.

hsize <- a %>%
  filter(Day > 49) %>%
  group_by(Run, I, C, R, W, gamma, Trajectory) %>%
  dplyr::summarize(H = h_entropy(MemoryEntropy) - H1,
            HippocampusDecrease = 100*(mean(MemoryEntropy)/H1 - 1),
            MeanTraumatic = mean(PTraumatic),
            NumTraumatic = sum(CTraumatic),
            MeanSimilarity = mean(ChunkSimilarity),
            Trajectory = )

## `summarise()` has grouped output by 'Run', 'I', 'C', 'R', 'W', 'gamma'. You can override using the `.groups` argument.

Visualizing Hippocampus Volume Decrease

For pretty labels, we define two new labeller functions:

ptev_val <- function(val) {
  paste("I =", val)
}

gamma_val <- function(val) {
  #expression(paste(gamma, "=", val))
  paste("gamma =", val)
}

hsize <- hsize %>%
  mutate(External = paste("I = ", I , ", C = ", C))

ggplot(filter(hsize, I != 1),
       aes(x=NumTraumatic, y=HippocampusDecrease, col = Trajectory)) +
  geom_point(alpha=0.1, size=0.1, pch=21) +
  #geom_density_2d() +
  facet_grid( W ~ gamma, labeller=label_both) +
  #scale_color_brewer(palette="Paired") +
  scale_color_viridis(option = "plasma", discrete = T, end = 0.8) +
  stat_summary(fun.data = mean_se, geom="point", size=1) +
  geom_smooth(method = "lm", formula = y ~ x, col="black", fullrange= F) +
  #geom_smooth(method = "lm", aes(fill=Condition), se=F) +
  ylab("Percentage Decrease in Hippocampus Volume") +
  xlab("Cumulative Number of Traumatic Memory Intrusions (Day 50-60)") +
  ggtitle("Symptom Severity And Hippocampal Volume Decrease") +
  theme_pander() +
#  theme(legend.position = "bottom") +
  theme(panel.background=element_rect(fill="NA", colour="black"))

By Trajectory

hsize_trajectory <- a %>%
  filter(Day > 49, I != 1) %>%
  group_by(Run, I, C, R, W, gamma, Trajectory) %>%
  dplyr::summarize(HippocampusDecrease = 100*(mean(MemoryEntropy)/H1 - 1),
                   PTraumatic = mean(PTraumatic),
                   CTraumatic = sum(CTraumatic),
                   MeanSimilarity = mean(ChunkSimilarity))

## `summarise()` has grouped output by 'Run', 'I', 'C', 'R', 'W', 'gamma'. You can override using the `.groups` argument.

And here is the corresponding plot:

ggplot(hsize_trajectory, aes(x=Trajectory, y = HippocampusDecrease, 
                             col=Trajectory, fill=Trajectory)) +
  facet_wrap(~I, labeller = label_bquote(
    rows = italic(I)[PTE] == .(I))) +
  scale_fill_d3() +
  scale_color_d3() +
  stat_summary(geom="point", fun.data="mean_sdl") +
  stat_summary(geom="errorbar", fun.data="mean_sdl", width=.25, 
               fun.args = list(mult=1)) +
  #geom_boxplot(width=.25, fill=NA) +
  geom_violin(alpha=0.4, col="NA", trim = T, width  = 1.5) +
  theme_pander() +
  ylab("% Hippocampus Volume Decrease") +
  theme(axis.text.x = element_text(angle = 90)) +
  theme(panel.background=element_rect(fill="NA", colour="black"),
        legend.position = "bottom")

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

In general, largest decreases of hippocampus occur in the Chronic trajectory, and the smaller decreases in the resilient ones. But is a larger decrease inverse;y associated with the proportion of resilient trajectories?

hsize_traj_prop <- hsize_trajectory %>%
  mutate(IsResilient = ifelse(Trajectory=="Resilient", 1, 0)) %>%
  group_by(I, C, R, W, gamma) %>%
  summarise(HippocampusDecrease = mean(HippocampusDecrease),
          PropResilient = mean(IsResilient),
          LogVolume = log(100 + HippocampusDecrease))

## `summarise()` has grouped output by 'I', 'C', 'R', 'W'. You can override using the `.groups` argument.

ggplot(hsize_traj_prop, aes(x=PropResilient, y=LogVolume, col=as.factor(I))) +
  geom_point(alpha=0.5, size=3, position = position_jitter(width=0.1)) +
  scale_color_viridis(option = "plasma", end=0.8, discrete=T) +
  geom_smooth(method = "lm", formula="y~x", col="black") +
  #scale_y_log10() +
  theme_pander()

ggplot(filter(hsize_trajectory, I != 1),
       aes(x=CTraumatic, y=HippocampusDecrease, col = Trajectory)) +
  geom_point(alpha = .2, size=0.5) +
  #geom_density_2d() +
  scale_color_d3() + 
  #scale_color_brewer(palette="Dark2") +
  #scale_color_viridis(discrete=T, option="inferno", end=0.8) +
  #stat_summary(fun.data = mean_se, geom="point", size=1) +
  #geom_smooth(method = "lm", formula = y ~ x, aes(col=as.factor(I)), fill="white",
  #            fullrange= T) + # theme_pander() +
  geom_smooth(method = "lm", aes(fill=Trajectory), 
              se=T, line_style="dashed") +
  scale_fill_d3() + 
  ylab("Percentage Decrease in Hippocampus Volume") +
  xlab("Total Number of Memory Intrusions (Day 50-60)") +
  ggtitle("Correlation Between Symptom Severity\nAnd Hippocampal Volume Decrease") +
  theme_pander() +
  theme(legend.position = "bottom") #+

## Warning: Ignoring unknown parameters: line_style

## `geom_smooth()` using formula 'y ~ x'

Functional Connectivity

k <- a %>%
  filter(! is.na(ChunkSimilarity),
         C == 0,
         Day > 50)

ggplot(filter(k, I > 1), 
       aes(x=Trajectory, y = ChunkSimilarity, 
           col=Trajectory, fill=Trajectory)) +
  stat_summary(fun.data=mean_se, geom="bar", alpha=0.5) +
  stat_summary(fun.data=mean_cl_boot, geom="errorbar", width=0.25, size=1) +
  facet_wrap(~I, ncol=3, 
             labeller  = label_bquote(
               rows = italic(I)[PTE] == .(I))) +
  ylab("PFC - Hippocampus Connectivity") +
  scale_color_d3() +
  scale_fill_d3() +
  theme_pander() +
  theme(axis.text.x = element_blank(),
        legend.position = "bottom")

Analysis of PTSD Model Simulations

Andrea Stocco

2/7/2021

Load and Transform the Data

Recovery Trajectories

Chronic Probability of Traumatic Intrusions

Finding a baseline to compare the model to existing literature

Factors affecting the distribution of trajectories

External factors: PTE Intensity and Context Similarity

Effects on Timecourses

Effects on Trajectories

Resilient Trajectory

Recovery Trajectory

Chronic Trajector

Delayed Trajectory

Internal Factors: Working Memory and Recollection Vividness

Effects on the Timecourse of Intrusive Memories

Effects on Trajectories

Resilient Trajectory

Chronic Trajectory

Recovery Trajectory

Delayed Trajectory

Comorbidty: Effects of Rumination

Effects on the Timecourse of Intrusive Memories

Effects on Trajectories

Neurobiology of PTSD

Hippocampus Size

Visualizing Hippocampus Volume Decrease

Functional Connectivity